Coulomb Drag Mechanisms in Graphene 
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Recent experiments revealed an anomalous Coulomb drag in graphene, which is distinct from the 
canonical mechanism relying on momentum transfer between graphene layers. A new mechanism, 
proposed to explain the observations, involves vertical energy transfer between layers, which couples 
to charge flow via lateral heat currents and thermopower. The two mechanisms depend on distinct 
physical effects, resulting in a starkly different behavior, in particular in drag magnitude and sign 
near charge neutrality. Under realistic conditions, energy transport dominates in a wide temperature 
range, giving rise to a universal value of drag which is essentially independent of the electron-electron 
interaction strength, providing a unique probe of energy transport on the nanoscale. The new 
mechanism explains strong enhancement of drag observed in classically weak magnetic fields. 



I. INTRODUCTION 

One of the most sensitive probes of interactions in 
low-dimensional systems is provided by measurements 
of Coulomb drag between adjacent electrically isolated 
electron systems, arising when current applied in one (ac- 
tive) layer induces voltage in the second (passive) layer. 
Coulomb drag was extensively studied in GaAs quantum 
wellsf^ where the observations were successfully inter- 
preted in terms of the momentum drag mechanisnPHSl 
(hereafter referred to as "P-mechanism" ) , in which inter- 
layer electron-electron scattering mediated by long-range 
Coulomb interaction transfers momentum from the ac- 
tive layer to the passive la yer. Recent measurements of 
Coulomb drag in graphene^ revealed behavior which 
clearly goes beyond the canonical P-mechanism. One 
novel aspect of graphene-based double-layer structures 
that makes them distinct from GaAs-based structures is 
that the typical separation of graphene layers can be as 
small as 1-2 nm, 8 which is much smaller than character- 
istic lengthscales such as the electron Fermi wavelength 
and the screening length. This defines a new strong- 
coupling regime wherein the interlayer and intralayer in- 
teractions are nearly equally strong. Fast momentum 
transfer between electron subsystems in the two layers 
and strong Coulomb drag have been predicted in this 
regime^HI 

Another new aspect of drag measurements in graphene 
is the unusual behavior observed in the double neutral- 
ity point (DNP) region. While away from DNP the ob- 
servations could be accounted for by models based on 
P-mechanism, it was found that drag exhibits a sharp 
peak in the DNP region which was not anticipated by P- 
drag models. These observations were explained by intro- 
ducing a new drag mechanism, which arises from inter- 
layer ener gy tr ansfer due to interlayer electron-electron 
scattering! 15 ! 16 ! xhe interlayer energy transfer can couple 
to charge currents, via lateral heat currents and ther- 
mopower, generating drag. Hereafter we will refer to 
this effect as "E-mechanism." It was predicted that in 
the DNP region the E-mechanism dominates over P- 



mechanism. Further, it was predicted that the si gn of 
drag at DNP is different for the two mechanisms! 15 * 16 ! 
These predictions provide a simple test for determining 
the leading mechanism in experiment. 

Somewhat paradoxically, while on a microscopic level 
the E and P mechanisms both arise from the same 
electron-electron interactions, the two contributions to 
drag are dominated by very different physical effects: 
momentum transfer between graphene layers via inter- 
layer electron-electron scattering vs. interlayer energy 
relaxation and long-range lateral energy transport cou- 
pled to charge flow. Notably, these effects arise on very 
different lengthscales. The length scales characteristic 
for the P-mechanism are on the order of the Fermi wave- 
length, which makes this mechanism essentially local. In 
contrast, the E-mechanism originates from lateral energy 
transport in the electronic system which is highly nonlo- 
cal. So far, the two mechanisms have been described by 
very different approaches, microscopic for P-mechanism 
and hydrodynamical for E-mechanism, which reflects the 
difference in the characteristic lengthscales. Here our aim 
is to develop a unified framework capable of describing 
these mechanisms on an equal footing. 

Given the difference in the characteristic lengthscales 
relevant for the E and P mechanisms, a suitable multi- 
scale framework is needed. This framework should also 
account for the peculiar dynamics of particles and holes 
near the Dirac point. First, electric currents carried by 
electrons and holes in the same d irecti on generate mo- 
mentum flow in opposite directions! 17 * 18 ! The same is true 
for energy flow due to particle and hole currents. In 
the presence of a magnetic field, the opposite sign of 
the Lorentz force on electrons and holes makes charge 
currents strongly coupled to neutral (energy) currents, 
resulting in a giant Nernst/Ettingshausen effectP^HlS] 

As we will show, this rich behavior is conveniently 
captured by a two-fluid model in which electrons and 
holes in the conduction and valence bands are described 
as separate subsystems coupled by mutual drag, origi- 
nating from both intralayer and interlayer carrier-carrier 
scattering. Here we will use a model analogous to that 



developed a while ago by Gantmakher and LevinsorPS 
to describe the peculiar magnetotransport phenomena 
in charge-compensated conductors, where magnetoresis- 
tance and Hall resistance exhibit anomalies at nearly 
equal electron and hole densities. A similar model was 
used to describe magnetotransport near CN in graphene 
in Ref!Ml Indeed, this model can successfully account 
for the strong influence of magnetic field on drag near 
DNP observed in Ref|Hl As we will see, P-mechanism 
and E-mechanism drag, for both B — and B^O, can 
be obtained from the same two-fluid model, allowing us 
to analyze these contributions to drag even-handedly. 

The two-fluid model predicts P and E contributions 
to drag with a very different dependence on density and 
magnetic field, which makes it easy to distinguish these 
contributions in experiment. This is illustrated in FigJTI 
where the E and P contributions are shown for several 
field values B. The E contribution features a large peak 
at DNP, whereas the P contribution is small near DNP 
and large away from DNP. The peak in the E contribu- 
tion is sharply enhanced by magnetic field, whereas the 
P contribution does not show strong field dependence. 
Overall, the density plots for P and E contributions look 
similar up to an overall sign reversal. Near DNP the 
E contribution is much larger than the P contribution, 
showing a distinct B dependence which closely resembles 
that reported in Ref 8 This behavior suggests 15 that the 
E contribution overwhelms the P contribution. We will 
analyze the relative magnitude of the E and P contri- 
butions below using the two-fluid model, and will show 
that under typical conditions the E contribution is indeed 
expected to be larger than the P contribution. 

The general reason for the smallness of the P contri- 
bution as compared to the E contribution can be under- 
stood as follows. The two-fluid model describes coupling 
between carriers of different types via the mutual drag 
coefficient 77, see Eq.((5]). The dependence of rj on the 
interaction strength ao = e 2 /hv can be obtained^ by 
matching the dependence of conductivity vs. ao found 
in Refs ll7|l"8l This gives a general relation of the form 
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where a is the RPA-screened interaction, N is spin/valley 
degeneracy and k is the dielectric constant for the 
substrate. Since N = 8 for the double layer system, 
the factor in the denominator is ~ 10 («o ~ 4 for a 
BN substrate). A ten- fold reduction of the bare value 
a sa 2.4 yields a small value of RPA-screened interac- 
tion, a « 0.25. The function F admits a power series 
expansion in a, arising from the solution of the quan- 
tum Boltzmann equation) 17 ^ 18 ! with the leading term be- 
ing a 2 ~ 0.06 (which corresponds to the two-particle 
Born scattering cross-section). This leads to a weak mu- 
tual drag, r) -C %. 

Crucially, the E contribution to drag remains unaf- 
fected by the small values of 77 so long as the interlayer 
thermalization occurs faster than the electron-lattice re- 
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FIG. 1: Coulomb drag originating from P and E mechanisms 
for different values of applied magnetic field, (a) Magnetodrag 
density dependence for P-mechanism, obtained from Eq. (171 
for B = 0.8 T, T = 200 K, 77 = 0.23ft. Note that drag resis- 
tance Pd,|| is weaker at DNP than away from DNP, and has 
positive sign at DNP. (b) Line traces along black dashed line 
in (a) for various values of B. (c) Magnetodrag density depen- 
dence for E-mechanism, obtained from Eq.(26l at B — 0.8 T, 
T = 200 K, p = 1. Note a large negative peak at DNP which 
is strongly enhanced by B field, a behavior distinct from that 
for P-mechanism. (d) Line traces along black dashed line in 
(c) for various values of B. Disorder broadening of DP of 
width A m 300 K was used in all plots (see Sec|V|. 



laxation. This is the case in graphene, since the electron- 
lattice cooling in this material is dominated by acoustic 
phonons, giving a slow cooling rate in a wide range of 
temperatures! 25 " 2 ' As a result, as discussed in more de- 
tail in Sec |IV[ the drag originating from E-mcchanism 
takes a "universal value" which shows little dependence 
on the interlayer scattering rate. To be more specific, as 
we will see in Sec |IV| the relative strength of the P and E 
contributions to magnetodrag at DNP can be estimated 
as 
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where a is the conductivity at charge neutrality and 
/3 ~ 1 is a factor describing temperature gradient buildup 
in respons e to energy flow in the system (see discussion 
IV below). Since a ~ Ae 2 /h whereas r\ is 
er than h, the factor in parenthesis is much 
smaller than unity. The smallness of the P contribu- 
tion, while being quite general, is not entirely universal. 
In particular, it does not hold far from DNP, where E- 
mechanism is small (see FigfTh. It also does not hold at 
elevated temperatures such that electron-lattice cooling 
length becomes small compared to system size, £ < W, 



leading to small values |3<1 which suppress the E con- 
tribution, see Eq.(25|. However, near DNP and at not 



too fast electron-lattice cooling, we expect E-mechanism 
to overwhelm P-mechanism. 

We also show that the anisotropy in energy transport 
due to geometrical constraints or contacts can have dra- 
matic effect on the E contribution, altering the density 
dependence. The difference in the behavior for the E 
and P contributions arises because these contributions 
are dominated by different physical effects. While the P 
contribution is controlled by electron-electron interaction 
strength, the E contribution is mostly controlled by the 
electron-lattice cooling length. In the regime when this 
length is comparable to system size, which is the case in 
few-micron-wide G/hBN/G devices, E-drag overwhelms 
P-drag over a wide temperature range. 

We will proceed as follows. First, we introduce the two- 
fluid model and employ it to describe the P-mechanism 
(Sees TTjTTT ). Next, we apply the two- fluid model to re- 
derive the E-mechanism, which allows us to compare the 
magnitudes and signs for the two mechanisms (Sec IV I. 



II. THE TWO-FLUID MODEL 

To describe transport near the Dirac point, it is cru- 
cial to account for the contributions of both electrons 
and holes. This can be done by employing the quan- 
tum kinetic equation approach of Refs ll7ll8l23l24l For 
a spatially uniform system, we have 



Qe(h) 
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dfe(h)(p) 

dp 



St[/ e (p),A(p)], (3) 



where f e (h)(p) is the distribution function for electrons 
and holes, and q e = —qu = e. The collision integral 
describes momentum relaxation due to two-particle col- 
lisions and scattering by disorder. 

The approach based on Eq.© is valid in the quasi- 
classical regime, when particle mean free paths are long 
compared to wavelength. This is true when the colli- 
sion rate is small compared to typical particle energy, 
which requires weak disorder 7 <C k&T, where 7 is de- 
fined in Eq. ( 35 1 , and weak effective fine structure con- 
stant a — e 2 Jhvo K "C I (k is the dielectric constant). 

The kinetic equatio n pi can be solved analytically 
in the limit of small a! 17 ! 1 ^. Rather than pursuing this 
route, we will adopt a two-fluid approximation used in 
RefsUn 24 which is particularly well suited for analyzing 
magnetotransport. In the two-fluid approach, transport 
coefficients can be obtained from the balance of the net 
momentum for different groups of carriers, electrons and 
holes, taken to be moving independently. We use a sim- 
ple ansatz for particle distribution function, 

/ e w(P)- e ( Ep -pa. w -l w )/ fcB r + 1 > ep = «»IPl. ( 4 ) 

where \x e = —fih are the chemical potentials of electrons 
and holes. The quantities a e and a/j, which have the di- 



mension of velocity, are introduced to describe a current- 
carrying state. This ansatz corresponds to a uniform 
motion of the electron and hole subsystems, such that 
the collision integral for the e-e and h-h processes van- 
ishes (as follows from the explicit form of the collision 
integral given in Reffl8|). Thus only the e-h collisions 
contribute to momentum relaxation, resulting in mutual 
drag between the e and h subsystems. 

In the following analysis we do not account for pos- 
sible temperature imbalance between electron and hole 
subsystems since fast e-e and e-h collisions quickly es- 
tablish thermal equilibrium locally in space. As we will 
see below, spatial temperature variation across the sys- 
tem becomes essential in the regime dominate d by energy 



1.15 



transport 10 . We will treat this regime in Sec IV 

Eq.© yields coupled equations for ensemble-averaged 
velocities and momenta of different groups of carriers, 
described by the distribution QpS 

- qi (E 4 + V, xB) = -^- r7VV(Vi - V*/), (5) 



where i, i' = 1,2,3,4 label the e and h subsystems in 
the two layers. The ensemble-averaged scattering times 
Tj, the carrier densities n,, and the electron- hole drag 
coefficient 77, describing collisions between electrons and 
holes, are specified below. The electric field Ej is the 
same for electrons and holes in one layer, but is in general 
different in different layers. 

The quantities Vj, Pi are proportional to each other, 
Pi — m,Vj. Here the "effective mass" is obtained by av- 
eraging over the distribution of carriers, as described in 
Sec Vj Eq.(33|. The integrals over p, evaluated numer- 
ically, give the effective mass as a function of T and /i. 
At charge neutrality, setting fi e th) = 0> we find 



9C(3) k B T 
2C(2) vl 



3.288- 



k B T 



(6) 



At high doping, /i 3> k^T, the effective mas s is given by 
the familiar expression, m — h/vq . In Sec III , we will 



use the approach outlined above to describe momentum 
drag. 

The two-fluid model can also be used to describe en- 
ergy transport. Indeed, particle flow is accompanied by 
heat flow, described by 



j q = TS e n e V e + TS h n h V h 



(7) 



where S e and Sh is the entropy per carrier for electrons 
and holes. Here the entropy and particle density can be 
related to the distribution function via 

s > = -— / ? f^[(i-/ i (p))Mi-/i(p)) (8) 

+/ i (p)ln/ i (p)] I ni = iJ^fi(p). 

In our analysis, we will need the value at charge neu- 
trality. Direct numerical integration in Eq.® gives 



S ~ 3.288fcB- In Sec IV we will connect j q to electric 
current, which will lead to a simple model for drag orig- 
inating from E-mechanism. 



III. MOMENTUM DRAG MECHANISM 

Here we will use the two-fluid model introduced in 
SeclITlto derive momentum drag. To facilitate the analy- 
sis of transport equations, it is convenient to switch from 
vector notation to a more concise complex-variable nota- 
tion. We will describe velocity, momentum and electric 
field by complex variables, 



V = V x + iV y , P = P x + iP y , E = E x +iE y . (9) 
In this notation, the transport equations take the form 



QiEi 



Ti 



iqiB\Vi = r)J2 n i'fe'-V^- ( 10 ) 



The solution of these equations can be written in a com- 
pact form by introducing the complex-valued quantities 



A; = 



iqtB + r]N ' 



N = 



i' = l..A 



(11) 



We can now express the velocity for each component as 
A; 



V, 



F + qi Ei 



F = rjy^mVi. 



(12) 
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FIG. 2: Hall drag originating from P and E mechanisms for 
the same parameters as in FigfT] (a) Hall drag density de- 
pendence for P-mechanism, obtained from Eq.(17l. (b) Line 



traces along black dashed line in (a) for various values of B. 
(c) Hall drag density dependence for E-mechanism , obtained 
from Eq.( |26| |. (d) Line traces along black dashed line in (c) for 
various values of B. The difference in sign for the two contri- 
butions makes them easy to identify experimentally. Disorder 
broadening of DP of width A ~ 300 K was used in all plots 
(see Sec(vf. 



Combining coupled equations for Vi, yields a relation 
(1 _ J2i' V^i')J2i n iVi = J2i QiMEi. For the quantity F, 
this gives 



F =?iE® A ^ A = E A * 



1-77A 



(13) 



Using this result, we evaluate the electric current in each 
layer by summing electron and hole contributions. For 
layer 1 we have 

jx = en le Vxe-en lh Vi h = eXx & {F+eEx)-eXxh{F-eEx)- 

(14) 
Combining this expression with a similar expression for 
layer 2, we can write the current-field relation for the two 
layers using a 2 x 2 matrix that couples variables in layer 
1 and layer 2: 



011 CT12 

021 °"22 



£l 

Eo 



(15) 



describe mutual drag between the layers (we note that 
012 = 021 )■ The real and imaginary parts of cri 2 describe 
the longitudinal and Hall drag. 

The longitudinal and transverse drag resistivity can be 
obtained by inverting the matrix, Eq.(15), giving 



Pd = Pd,|| + JPd.Hall 



012 



011022 — 012021 



(17) 



The quantities p d 11 and pd.Haii give the magnetodrag and 
Hall drag shown in Figs Tj2 This quantity features an 
interesting dependence on carrier density and magnetic 
field. We will analyze the limit of small 77 (weak interac- 
tions). In this case, we have A^ rj «j/( — — iq%B). This 
gives a xx ~ e 2 (X le + Xxh), 022 ~ e 2 (A 2e + X 2 h), which 



describes conductivity of decoupled layers, 
terms are small in r\: 



The cross- 



012 = 021 ~ e 2 ?7(Ai e - Ai /l )(A 2e - X 2h ) (18) 



Here oxx = e 2 [(A le - X lh )fx + X le + X lh ], 012 = 021 = This S ives the dra § resistance 



e 2 (Ai e - Xxh)f2, 0-22 = e 2 [(A 2e - X 2h )f 2 + A 2e + A 



21, \ 



h 



rj{X le - Xi h ) 
1-rjA : 



, ??(A 2e -A 2ft ) 
h ~ l-r,A ■ (16) 



Pd 



->l 



(Ai e — Aih)(A 2e — A 2 h) 
e 2 (A ie + Xxh)(X 2e + X 2h ) 



(19) 



Here the quantities oxx an d ct 22 describe the conductivity 
of layers 1 and 2, whereas the quantities 0x1 and a 2 x 



For B — 0, this quantity vanishes at the double neutrality 
point (DNP), nie = nxhi n 2e = n 2 h- Drag is negative for 
equal-polarity doping and positive for opposite-polarity 



doping, reproducing behavior well-known for momentum 
drag. 

For nonzero B, the Hall drag and magnetodrag can 
be obtained b y ex panding ImAi(-B) = Xi(0)(qiTi/mi)B + 
0(B 3 ) in Eq.(19|. This gives Hall drag vanishes exactly 
at DNP but is nonzero near DNP. For electron and hole 
densities near DNP, such that A e ~ A/ l; we find 



PcLHall 



)]T 



B 



(Ai e — Alfc) (A 2e — A 2 /j) 



0(B 3 



em \(Xie + Xih) (A 2e + A 2 /0 

This expression vanishes on the line n\ = — n 2 corre- 
sponding to doping of opposite polarity in the two layers. 
In contrast to Hall drag, magnetodrag is nonzero at 
DNP. From Eq.(19l we obtain a finite magnetodrag of a 
positive sign: 



Pd,\\ = V- 



t 2 B 2 



0{B 4 



(21) 



Here the quantities r and m are evaluated at charge neu- 
trality, n e = n/j, in each layer. Interestingly, the mag- 
netodrag sign comes out opposite to the sign predicted 
by the energy transport modelP^ The magnetodrag sign 
therefore provides a clear signature which discriminates 
between the E and P mechanisms in experiments. 

The density dependence of magnetodrag and Hall drag 
predicted from P-mechanism, which is shown in Figs Il2 



is obtained by numerically evaluating the expression in 
Eq.([T7| [see Secfvlfor detailed discussion of the modeling 
procedure]. In agreement with the above analysis, p^ 
in Fig. [Ik,b is positive at DNP, increasing quadratically 
with B field. Also, pa.Haii m Fig. [2k, b increases linearly 
with B field vanishing along n\ — — n 2 as expected. 



IV. ENERGY-DRIVEN DRAG MECHANISM 

To compare the above results to those obtained from 
the energy-driven drag mechanism introduced in Ref ll5l 
it is useful to re-derive the latter using the two-fluid 
model. First, we need to relate the quantity j q [Eq.([7])] 
and electric current. It will be instructive to ignore 
the mutual drag effect discussed above, setting r\ = 0, 
and restore finite r\ later. In this case, continuing to 
use complex variables for velocities and fields, we find 
V e = (X e q e /n e )E, V h = (X h q h /n h )E. Applying this to 
Eq. ([T]) yields a relation 



jq = Qj, Q = 



T (S e X e q e + S h X h q h ) 



A,. 



X h q\ 



(22) 



where S e nA can be evaluated using the expression in Eq. 
[8| This relation is particularly useful since the effect of 
the Lorentz force is fully accounted for via A^. Specif- 
ically, the imaginary part of Q generates relative phase 
between complex variables j q and j, which describes the 
Nernst/Ettingshausen effect. 

Energy transport, described by Eqs.(|7|),(|22"l), induces 
a buildup of temperature gradient across the system. 



For two layers in close proximity, fast heat transfer be- 
tween layers due to interlayer electron scattering leads 
to a temperature profile which is essentially identical in 
both layers^. The temperature gradients can drive a 
local thermopower via 



Q 

E=^VT 



(23) 



where the quantity Q is given by the ratio of the heat 
current and electric current for the layer in question. As 
discussed in detail in RefQ~5] this relation follows from 



Onsager reciprocity combined with Eq.(22). 



The temperature gradient can be found from balancing 
the heat flux due to thermal conductivity against the net 
heat flux in the two layers. While the details of the anal- 
ysis somewhat depend on sample geometry (see Ref ll5l 
and discussion below), we can write the general relation 
for temperature gradient in a general form as 



( Kl + K 2 )vr = /3(j 1 , q + j 2 , q ) 



(24) 



where k± and k 2 are thermal conductivities of the layers, 
and the factor < j3 < 1 describes the "active part" 
of the heat that is not lost to contacts and the crystal 
lattice. 

The value of (3 depends on the relation between 
electron-lattice cooling length and system dimensions. 
For example, for the case when cooling to the lattice 
dominates heat loss, Ref. 15 finds the value 



P 



1 



1 



tanh c 



W 
2£' 



(25) 



where W is system size and £ is the electron-lattice cool- 
ing length. In this case, (3 —>■ 1 when £ ^> W and j3 — > 
when £ <C W. The cooling length, £, in graphene can be 
as large as several microns for a wide range of tempera- 
tures up to room temperatur er 5 " 27 [ As a result, since in 
typical devices W is a few microns or smaller, the factor 
/3 can be close to unity j3 ~ 1. However, at tempera- 
tures exceeding a few hundred kclvin the cooling length 
£ shortens rapidly owing to cooling pathway mediated 
by optical phonons. At such high temperatures Eq.(25) 
predicts vanishingly small (3. In this case, temperature 
gradients in the electron system do not build up. Here 
we concentrate on the technologically important temper- 
atures T < 300 K that result in appreciable j3 values. 

The quantity /3 can in principle be obtained by solv- 
ing the heat transport equations, as in Ref. 15. Here, 
however, we will focus on general behavior that can be 
understood directly from Eq.(24) without specifying /3 
value. In the isotropic case, where heat flow is not influ- 
enced by device geometry or contact placement, /3 is a 
real number (in the anisotropic case (3 becomes a tensor, 
see below). Treating f3 as a number, the drag voltage can 
be found by combining the above relations to obtain 



P12 =13 
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K2j 



(26) 



This expression is symmetric under interchanging lay- 
ers, 1 -H- 2, and features interesting dependence on car- 
rier density shown in Fig{TJ:,d and Fig(2] c,d. Notably, 
the signs of magnetodrag and Hall drag obtained from 
E-mechanism are opposite to those obtained from P- 
mechanism. This provides a convenient way to differen- 
tiate between the two mechanisms in experiment. Here 
we used device parameters indicated in the caption [see 
Sec jV] for the discussion of the modeling procedure]. 

The behavior is particularly interesting at charge neu- 
trality where, since n e = Uh, S e = Sh, the two terms 
contributing to the heat current j q are of equal magni- 
tude. In this case, because A e = A^ at zero magnetic field 
B, drag resistivity vanishes at B = 0. Importantly, since 



A e = A^ at finite B, the quantities Qi that enter Eq.(26) 



are purely imaginary. As a result, drag resistivity has a 
negative sign for nonzero B. We can obtain magnetodrag 
by expanding in small B, which gives 
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FIG. 3: Long-range energy flow which governs E-mechanism 
leads to anisotropy in the quantity f3 which relates heat flow 
and temperature gradient, Eq.(30|. The effect of system 



anisotropy on drag originating from E-mechanism is illus- 
trated for /3 given in Eq.(31l. Density dependence of Pd,|| 
and pd.Haii, Eq.fMl, is shown for B = 0.8 T, T = 200 K and 



other parameter values identical to those used in Figs|l|2 
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where r, m and k are evaluated at charge neutrality. For 
an estimate, we will relate thermal conductivity to elec- 
trical conductivity using the Wiedemann-Franz relation, 

k = „ f a. This relation is valid for degenerate Fermi 
system, however we expect it to be also approximately 
valid near charge neutrality. This gives 
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Comparing to the answer for P-mechanism, we find the 
ratio of the contributions due to momentum and energy 
mechanisms 



(P) 
« (E) 



2k 2 



ay 



zps 2 



We can estimate entropy per carrier at DP by evaluat- 
ing the integral over energy in Eq.|8]). Using the value 
S ~ 3.288/cb quoted above, we arrive at Eq.pl). Given 
the conductivity value at charge neutrality, a ~ 4e 2 //i, 
and taking into account that the mutual drag coefficient 
r\ is s mall when the fine structure constant a = e 2 /Tivk is 
smalP^21, we conclude that the ratio in Eq. Q is smaller 
than unity. This indicates that under very general condi- 
tions the E contribution overwhelms the P contribution 
in the DNP region. 

Since E-mechanism depends on long-range energy 
transport, it is particularly sensitive to sample geometry 
and the arrangement of contacts. As discussed in Ref. 
[T5T the source and drain contacts can act as ideal heat 
sinks that suppress temperature gradients along current 
flow. In contrast, for a Hall bar with "non-invasive" volt- 
age probes, temperature gradient in the direction trans- 
verse to current flow will remain essentially unaffected by 
the voltage probes. This anisotropy can be easily imple- 



mented by making fl in Eq.(24| a tensor, 



$■ 



/ Pxx Pxy 

Pyx Pyy 



(30) 



so that heat current in one direction generates higher 
temperature gradient than heat current in another di- 
rection. For example, the Hall bar geometry analyzed 
in Ref. [15] is described by a tensor with eigenvalues ap- 
proximately equal and 1 for the directions along and 
transverse to the Hall bar, respectively. Choosing the x 
axis along the bar, we obtain 



& 




1 



(31) 



Plugging this matrix in Eq.(24) and proceeding as above, 



(29) we arrive at an expression for drag resistivity 



P12 



iQ 2 lm(Qi) 
T(ki 
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(32) 



The dependence on density can be obtained by plugging 
Q from Eq.(22| in EqJ32|). This gives pd,|| and /3d,Haii 
maps shown in Fig. [3^75. Note the absence of layer 
symmetry 1 -H- 2 which was present in the expression 
(1261) and manifest in drag maps in Figs Tl2 The density 



dependence shown in FigJ3]is essentially identical to that 
found in ReflT51 



V. MODELLING PROCEDURE 

Here we comment on the quantities that enter the two- 
fluid description, and discuss the sensitivity of the results 
to the simplifying assumptions made in the model. 

In the two-fluid model we describe the momentum- 
velocity relation for each component as P = m^V, where 
rrii is an "effective mass." An explicit expression for m t 



as a function of T, u can be found by expanding the dis- 
tribution functions Q to lowest non-vanishing order in 

a e(h)'- 



i J^pp^^Mp) _ 1 Jd 2 p P l 9l (p) 



VofcPpfVaJiip) vofdipfg^p)' 



(33) 



where &(p) = /,(p)(l - /*(p)). 

The times r^ for disorder scattering and carrier densi- 
ties rii in Eq.([5]) are expressed through the distribution 
function p} with a; = 0: 



d 2 P MP) 
(^) 2 n(e p y 



m = 4 



d 2 P 



/i(p), (34) 



where the factor of four accounts for spin-valley degen- 
eracy in each layer. We pick a model for the transport 
scattering time r(e) to account for the experimentally 
observed linear dependence of conductivity vs. doping, 
o = jJL*\n\, where //* is the mobility away from the DP. 
This is the case for Coulomb impurities or strong point- 
like defects, such as adatoms or vacancies^. In both 
cases the scattering time has an approximately linear de- 
pendence on particle energy, 



00 
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1 = v 



y/eh/fj,* 



(35) 



where the disorder strength parameter 7 is expressed 
through mobility. The value /i* = 6 ■ 10 4 cm 2 /V • s mea- 
sured in graphene on BN- 9; yields 7 « 120 K. Similar val- 
ues for 7 are obtained from the DP width extracted from 
the resistivity density dependence,™ An « 5 • 10 10 cm -2 . 
In doing simulations, we found it convenient to use a 
different, simplified model for transport scattering which 
does not involve integration over particle distribution, 
yet yields results similar to those obtained from a more 
microscopic model, Eq.(34). We model the scattering 



time in the full range of doping densities as 



= (rriiV 2 + A) 



A 2 



(36) 



where m, depends on temperature and density via 
Eq. ( 33 1 , and the parameter A describes the smearing of 
DP due to disorder. This model accounts for the ex- 
perimentally observed linear dependence of conductiv- 
ity vs. doping. In the simulation we used the value 
A = 300 K which translates into DP width of about 
An « 10 11 cm -2 , consistent with the above estimates. 

Additionally, we found it convenient to account for dis- 
order broadening of DP by using an effective tempera- 
ture T ff = T + A in the evaluation of the effective mass 
in Eq.(33) and the entropy per particle in Eq.®. This 



simple procedure captures the essential characteristics of 
DP broadening since smearing of the density of states by 
temperature and disorder occur in a similar fashion. 
Using the parameters A = 300 K, T = 200 K and 



the model for scattering time Tj in Eq.(36l, we plot p^u 
and pd.Haii for P-mechanism [Eq.(17l] and E-mechanism 
[Eq.(26l] in Figs lp respectively. As discussed above, 



this gives density dependence of drag that differs in sign 
for the two mechanisms. Additionally, we find that the E- 
mechanism magnitude exceeds that of P-mechanism for 
the region near DNP. This agrees with the small ratio for 
P-mechanism vs. E-mechanism derived for small n. 



VI. SUMMARY 

In summary, we argue that drag in graphene near 
charge neutrality is dominated by energy transport ef- 
fects (E-mcchanism) arising due to fast interlayer en- 
ergy relaxation that couples to lateral energy flow and, 
via thermopower, drives electric current. We devel- 
oped a two-fluid framework which accounts both for the 
E-mechanism and for the standard momentum-transfer 
drag (P-mechanism), capturing the essential features of 
the two mechanisms. This unified approach is particu- 
larly instructive, not only because it produces both P- 
mechanism and E-mcchanism, but also because it allows 
an unbiased way of comparing the magnitudes of the two 
mechanisms. Strikingly, the P and E mechanisms yield 
opposite sign for both magnetodrag and Hall drag re- 
sistivities. Along with a strong peak in magnetodrag at 
DNP originating from E-mcchanism, this provides a clear 
way to experimentally distinguish the two mechanisms. 

We show that the magnitude of drag originating from 
the two mechanisms are dominated by very different ef- 
fects. The P mechanism is mostly controlled by the inter- 
layer electron-electron interaction, becoming weak when 
this interaction is screened. In contrast, the E mechanism 
is controlled by long-range energy transport, yielding a 
"universal value" for drag which is essentially indepen- 
dent on the interaction strength. Estimating typical de- 
vice parameters, in particular the electron-lattice cooling 
length, we conclude that drag at DNP is dominated by 
the E mechanism. Drag measurements in graphene can 
therefore be used as a unique probe of energy transport 
on the nanoscale. 

We acknowledge useful discussions with A. K. Geim, 
P. Jarillo-Herrero, L. A. Ponomarenko, and financial sup- 
port from the NSS program, Singapore (JS). 

Upon completion of this manuscript, we became aware 
of a related work by Titov et al.™ 
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